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ABSTRACT 

The flux of ISM neutral atoms surrounding stars and their environment affects the 
motion of dust particles in debris disks, causing a significant dynamical evolution. 
Large values of eccentricity and inclination can be excited and strong correlations 
settle in among the orbital angles. This dynamical behaviour, in particular for bound 
dust grains, can potentially cause significant asymmetries in dusty disks around solar 
type stars which might be detected by observations. However, the amount of orbital 
changes due to this non-gravitational perturbation is strongly limited by the coUisional 
lifetime of dust particles. We show that for large values of the disk's optical depth the 
influence of ISM flow on the disk shape is almost negligible because the grains are col- 
lisionally destroyed before they can accumulate enough orbital changes due to the ISM 
perturbations. On the other hand, for values smaller than 10^"^, peculiar asymmetric 
patterns appear in the density profile of the disk when we consider 1-10 mum grains, 
just above the blow-out threshold. The extent and relevance of these asymmetries 
grow for lower values of the optical depth. An additional sink mechanism, which may 
prevent the formation of large clumps and warping in the disks is related to the fast 
inward migration due to the drag component of the forces. When a significant eccen- 
tricity is pumped up by the ISM perturbations, the drag forces (Poynting-Robertson 
and in particular ISM drag) drive the disk particles on fast migrating tracks leading 
them into the star on a short timescale. It is then expected that disks with small 
optical depth expand inside the parent body ring all the way towards the star while 
disks with large optical depth would not significantly extend inside. 
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1 INTRODUCTION 

The observable part of debris disks are small (^J 1 mm) 
dusty or icy grains, coUisionally produced from larger, un- 
detectable parent bodies. In addition to the gravitational 
pull of the star, these grains are also affected by several 
forces such as stellar radiation pressure, Poynting-Robertson 
(PR) drag and the possible gravitational influence of large 
bodies in the neighborhood. As has been shown in numer- 
ous numerical studies, the combined effect of these differ- 
ent forces can le ad to comple x spatial structures in re- 
solved disks (e.g. I WvattI [20081 ') . A less investigated addi- 
tional force that could have an influence on grain dynamics 
is the drag due to particles from the surrounding interstellar 
m edium (ISM). The effect of ISM has first been addressed 
by I Artvmowicz fc ClampinI (|l997l ). who studied the level of 
disk erosion due to sandblasting by ISM dust grains. They 
concluded that, at least around massive stars, this effect was 
negligible because small ISM grains felt a strong repulsive 



radiat i on force. More recen tly Scherei ll2000ll, Debes et al.l 
ll2009l) [MTness et al.l l|2009l ). lBelvaev fc Rafikovl l|2010l ') aiid 
'Pastoil ( 201ll ) considered instead the effect of ISM neutral 
atoms on disk grains. This flux of neutral atoms acts in- 
deed similarly to the solar wind or radiation pressure from 
a physical point of view but, being monodirectional, can 
significantly perturb the trajectories of the grains, and po- 
tentiall y induce asymmetric struc tures in the d isk. In par- 
ticular iManess etafl (|2009l ') and iDebes et al] (|2009l 'l sug- 
gest that the ISM flux can explain the unusual morphol- 
ogy of some debri s disks like HD61005 and HD32997. In 
their model Debes et al.l (|2009l ) consider dust particles close 
to the blow-out size for the star and compute the trajecto- 
ries of perturbed grains over a timescale of 5000 yrs. The 
majority of their grains are strongly perturbed and end up 
quickly on hyper bolic orbits. A similar scenario is outlined 
by Maness et al. I |2009) where they concentrate on small 
grains (0.1 ^m) whose lifetime before ejection is of the order 
of a few 10^ years. The morphology changes they observe is 
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mostly due to the fast transfer of grains from low eccentricity 
orbits into hyperbolic trajectories. 

In this paper we take these studies a step further and we 
concentrate on the effects of the ISM neutral flow on bound 
Keplerian orbits of dust particles in debris disks around so- 
lar type stars. The orbital changes on relatively large dust 
grains (we model grains with radius ranging from 1 to 10 
/xm) is slow and it occurs on timescales of the order of 10® 
yrs. This timescale is m uch longer than that considered in 
iDebes et al] l|2009f) and iManess et ahl l|2009l ') and it is re- 
lated to the larger size of the grains. In addition, our grains 
are bound to the star and the orbital evolution induced 
by the ISM flow lead to eccentric and possibly asymmet- 
ric disks, characterized by density clumps, only if the col- 
lisional lifetime of the grains allows it. Thus, a crucial pa- 
rameter that controls the efficiency of the ISM flow in shap- 
ing debris disks is their geometrical optical depth, on which 
coUisional lifetimes directly depend. The coUisional cascade 
that is steadily eroding, by cratering and fragmentation, all 
solid bodies in a debris disk, does indeed reduce the life- 
time of dust grains, thus limiting the amount of time their 
trajectories can be perturbed by non-gravitational forces 
like radiation pressure, PR drag and interaction with ISM 
l|Wvattll2005l ). To correctly evaluate the impact of ISM on 
the density profile of a dust disk we thus need to account for 
the limited lifetime of individual particles and the amount of 
orbital changes that they can accumulate during that time. 

The rol e of collisions had bee n discarded bv lDebes et al.l 

l|2009l) andlM aness et al.l (|2009l ) because their dust grains 
are pushed on hyperbolic trajectories by the ISM flow on a 
very short timescale (a few 10'^ years) , comparable to the col- 
lisional lifetime of the particles they handle. In our scenario, 
grains remain on bound orbits but a simple numerical inte- 
gration of dust trajectories over a long timespan, performed 
to predict the overall disk density distribution, can lead to 
misleading result which are correct only for disks with very 
low optical depth. A numerical model intended to evaluate 
the impact of ISM on the morphology of dense disks must 
include an estimate of the lifetime of each individual parti- 
cle that will be part of the disk only for an interval of time 
no longer than its coUisional lifetime after which it will be 
replaced by a new one. If the integration time of the whole 
particle ensamble is a few times longer than the lifetime of 
the considered dust particles, we will obtain a stationary re- 
laxed population of grains which will describe in a reliable 
way the spatial distribution of the debris disk. We must also 
account that dust of different sizes would have different life- 
times so the presence of density structures would depend on 
the size of the dust that is imaged. 

The effect of ISM neutral atom flux on the dust grain 
trajectories is twofold. On one side it forces a coupled evo- 
lution of eccentricity and pericenter longitude that drives 
the particles on aligned eccentric orbits. On the other side, 
if the ISM flow is inclined with respect to the parent body 
orbital plane a significant inclination can be pumped up. In 
addition, a fast inward migration is forced by the PR drag, 
which is strong for highly eccentric orbit, and by the ISM 
drag component of the force related to difference between 
the monodirectional velocity of the ISM neutral atoms and 
the orbital velocity of the dust particles. When they mi- 
grate close to the star, the grains either sublimate or impact 
potential planets or the star itself going lost. This is an addi- 



tional sink mechanism, purely dynamical, which contributes 
to the erosion of the dust population. 

In this paper we intend to address these issues con- 
centrating on the balancing between the orbital evolution 
of dust grains in dusty disks under the effect of ISM neu- 
tral atoms perturbations and coUisional lifetime. Our goal 
is to outline the features of the steady state population of 
dust produced by the combination of these effects. We will 
explore the amount of asymmetry in the debris disk cre- 
ated by the ISM perturbations and how this depends on the 
optical depth of the disk. We concentrate in particular on 
dust grains of 1 um size since these particles, which are just 
above the cut-off size imposed by radiation pressure(~ 0.6 
^m), are those that dominate the geometrical cross section 
of the system, and thus the flux comi ng from the disk in the 
visible and near-IR (see for example iThebault fc Aueereaul 
(2007)). We will also give a hint to the evolution of 10 um 
size particles to test how the density distribution of the disks 
depends on the observed particle size. 



2 THE NUMERICAL MODEL 

Our numerical integration scheme computes the trajecto- 
ries of a l arge number of dust particles of a given size using 
RADAU (jEverhartlll985h . It is a very precise algorithm that 
accurately handles highly eccentric orbits as those produced 
by ISM perturbations. It is an implicit Runge-Kutta in- 
tegrator of 15th-order which proceeds by sequences within 
which the substeps are taken at Gauss-Radau spacings. 
High orders of accuracy are obtained with relatively few 
function evaluations. In addition, after the flrst sequence, 
the information from previous sequences is used to improve 
the accuracy and the integrator itself chooses the next se- 
quence size depending on the estimated variability of the 
'force' term in the differential equations. Non-gravitational 
fo rces like radiation pressur e and PR drag are included as 
m 

Im 

arzari fc Vanzanil l| 19941 ) and, in addition, we model the 
effects of the ISM neutral gas atoms adopting the s ame ap- 
proach used to calculate the solar wind pressure (|Schereij 
(2000)). The acceleration of a dust particle due to the im- 
pact with fast neutral atoms of the ISM wind is computed 
as: 

f = -~CD{nHmH)As I V - vh I (v - vh) (1) 

where v is the orbital velocity of the dust grain and 
Vh is the velocity of the ISM neutral atoms, mn and nn 
are the mass and density of neutral hydrogen atoms, respec- 
tively. As = 'KT^ is the area of the dust particle of radius r 
while Cd is a drag coefficient whose value is about 2.5. This 
functional dependence is robust only in the limit where the 
gas of the ISM has a temperature lower than about 5000 K. 
Under this condition, the thermal speed is lower than vh 
and the momentum change is proportional to the {v — vhY ■ 
For hot (T > 5000^!") ISM cloud s the momentum change 
would be proportional to {v — vh aev fc Rafikovllioiol : 

iDraine fc Salpeterlll979h . 

For the coUisional lifetime of each particle, we adopt the 
usual expression 

T 
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where r is the local geometrical vertical optical depth of 
the disk and T the orbital period of massless particles at this 
location (e.g. lArtvmowicz fc Cla mpin 1997) . Although de- 
tailed numerical studies fe.g.. lThebault fc Augereau (,2007 )) 
show that real collision rates may strongly depart from this 
simplified expression, we adopt it here for its simplicity and 
because it is an easy way to tune in the coUisional effects. 

Particle size is a crucial parameter since it determines 
the response to radiation pressure parameterized by j3, the 
ratio of radiation pressure to stellar gravity. Given the steep 
size distribution for coUisional cascades and the limited num- 
ber of test particles in a deterministic code, it is impossible 
to consider a wide size range. We thus restrict our study to 
several independent simulations each with single size parti- 
cles. These different runs can then be combined, when prop- 
erly weighted, to give a first-order estimate of "real" disks. 
We consider a scenario where the central body is a solar type 
star and the dust grains of the debris disk are the by-product 
of collisions occurring in a ring of parent bodies (planetesi- 
mals, asteroids, comets) on circular orbits. Whenever a dust 
particle is produced, its orbital elements are computed from 
the position and velocity vectors of the parent body account- 
ing for the reduced gravity 1-/3 due to radiation pressure. 
The average coUisional lifetime tc is derived for the particle 
and, at starting, a random value of 'age' ta is given to each 
grain. As the integration advances, ta is updated with the 
timestep At. When ta becomes larger than the coUisional 
lifetime tc, the particle is eliminated from the sample and a 
new one is generated with new orbital elements and ta = 0. 
When the timespan of the integration is a few times longer 
than tc we have a steady state population of dust parti- 
cles where the coUisional lifetime determines for how long a 
grain is subject to non-gravitational perturbations. In addi- 
tion, when a particle is injected on a hyperbolic orbit or it 
has migrated far inside the parent body ring it is discarded 
and a new one is generated. As a consequence, the steady 
state population can be due to a balance between the differ- 
ent mechanisms of dust elimination i.e. collisions, ejection 
out of the system or inwards migration. We assume in our 
model that the parent bodies are large enough not to be sig- 
nificantly affected by ISM. This is a reasonable assumption 
since the forces acting on the particles strongly depend on 
their size. Larger parent bodies would be affected by ISM 
only on a much longer timescale. 



3 INITIAL SETUP 

We consider as a test case for our modeling a disk produced 
by a ring of parent bodies evenly distributed in between 50- 
70 AU on circular orbits. Since the size of the ring is small, 
assuming a uniform distribution is a good approximation. 
We adopt different values for the optical depth of the disk 
ranging from r ^ 10""^, typical of dense collision-dominated 
disks like that around Beta Pictoris, to r ~ 10~^ for almost 
coUisionless disks similar to the Kuiper Belt where trans- 
port mechanisms like PR drag, stellar or ISM winds can 
dominate. The ISM gas flow is approximated as a monodi- 
rectional flux moving along the x-axis with a velocity of 20 
km/s and a concentration of neutral hydrogen atoms equal 
to uh = 0.1 cm~^. These are typical values of the inter- 
stellar cloud surrounding our solar system. The local speed 



of the sun with respect to the ISM is around 25 ± 2 km/s 
but some deceleration of the interstellar neutral hydrogen 
flow across the heliospheric interface is expected reducing 
the speed to about 22 km/s (LaUement 2Q0^. The density 
we adopt in the numerical m odel is intermediate between 
the value given bvlFahJ (Il996l) e qual tonfr = 0.05 cm and 
that given in lFrisch et all (|l999l ) of uh = 0.22 cm"^. 



4 THE COPLANAR CASE: {E,txj) EVOLUTION 

We first assume that the interstellar gas flow lies on the 
plane of the parent body ring. This may be considered a 
test bench where to study the evolution of the (e, zj) or- 
bital elements. In effect, it may not be a realistic scenario 
since the fiow, impinging on one side of the disk, may not 
reach with equally intensity the other side being eventually 
absorbed. This does not occur when there is some inclina- 
tion between the dust orbits and the ISM flow. However, 
it is useful to explore the density distribution in the planar 
case since the structures and dumpings that develop in the 

(x y) plane are retrieved when a low inclination between 

the flux and the initial dust orbits is assumed. In Fig[l]the 
position of 1 jim size particles are given after an evolution 
lasting about 3 times the average coUisional lifetime for the 
case with r = 1 x 10"^ (At = 0.5 Myr) and r = 1 x 10"* 
(At = 5 Myr, about twice Tstark defined by formula (7)). In 
the case with r = 1 x 10"^ the coUisional timespan is around 
100 Myr, however the fast inward evolution due to the large 
eccentricity, achieved by the dust grains because of the ISM 
perturbations, significantly limits the lifetime of the particle 
within the disk. The most relevant sink mechanism in this 
scenario is related to the fast inward migration of grains 
driven by both the PR drag and the ISM drag components 
of the force. 

The density plots are derived by computing the num- 
ber of particles populating the local spatial area (equal size 
squared bins) at the end of the numerical simulation. This 
number is then normalized to the total number of particles 
which is kept constant and equal to 4 x 10* in all our simu- 
lations. 

The first plot on the top left in Fig[l] shows how the de- 
bris disk appears when the ISM force is neglected. It is axis- 
symmetric with a density peak at the internal ring corre- 
sponding to the location of the parent bodies. This is where 
all the pericenters of the particles reside and the higher den- 
sity is due both to the geometrical configuration (the orbit 
is tangent to the circle) and to the creation of new grains 
when the older ones are destroyed. A second peak in the 
density plot appears at the outer ring where the apocenters 
of the particles are located and the grains move at the lowest 
speed. The pericenter longitudes are randomly distributed in 
between and 360° when the steady state is reached. The 
initial eccentricities of our sample of particles are encom- 
passed between 0.42-0.45 while the initial semimajor axes 
range from 90 to 120 AU. This orbital distribution is due to 
the value of /3 adopted at the beginning of the simulation 
that reduces the mass of the central star injecting the grains 
into eccentric orbits. 

When the ISM effects are included, asymmetries de- 
velop in the disk due to the eccentricity and pericenter evo- 
lution forced by the ISM perturbation. These asymmetries 
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Figure 1. 2-D density distribution of 1 /jm size dust particles in debris disks generated by a ring of parent bodies moving on circular 
orbits in between 50-70 AU. The density is computed as number of particles populating a squared region at a given timestep divided by 
the total number of particles in the model. The top-left plot shows how the disk appears when the ISM flow on the grains is neglected. 
The top-right plot illustrates the density distribution when ISM is acting on the particles and the optical depth is r = 1 X 10"'^. The 
two bottom plots show the density distribution when t = 1 X 10~* and r = 1 X 10~®, respectively. 



are barely detectable when r = 10""^ (Figdl upper right 
plot) but they become more pronounced for lower values of 
r (Fig[T] bottom plots) since the grains have more time to 
accumulate the perturbative effects. When r = 10~^ the two 
peaks at the particles pericenters and apocenters are still vis- 
ible but the disk begins to develop an eccentric shape. The 
alteration in the disk shape is more noteworthy for r — 10~® 
when the disk appears as the superposition of two ellipses. 

This feature can be explained by inspecting the orbital 



behaviour of the grains under the effect of the ISM neutral 
gas. In Fig[2]we show the evolution with time of the grain 
orbits in the (e, m) plane. Particles evolve towards high ec- 
centricity values (~ 1) while, at the same time, the pericen- 
ter values tend to 270°. When large values of eccentricity 
are achieved, the particles drift quickly inside because of 
the drag component of the forces. 

This behaviour can be interpreted on the basis of 
the analytical solution of the classical Stark problem 
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l|Belvaev fc Rafikovllioiol : |Pastorll20lil ') where the Keplerian 
motion of the particle is perturbed by an external force con- 
stant in magnitude and direction. In effect, the dust particle 
dynamics is different from that predicted by the classical 
Stark problem since it includes PR drag which is a tangen- 
tial force. In addition, the ISM force depends on the differ- 
ence between the particle v and the neutral atom velocity 
vh so that the force changes periodically depending on the 
position of the particle in its orbit. In our scenario, v is 
approximately 5-10% of vh but it can become as large as 
30% for highly eccentric orbits once close to perihelion. This 
has a significant effect on the semimajor axis of the particle 
orbit causing a fast inward migration. However, the analyt- 
ical equations that describe the evolution of a body in the 
general Stark problem gives a very accurate description of 
the evolution of eccentricity, pericenter longitude, inclina- 
tion and nodal long i tude. The secular solution, obtained by 
iBelvaev fc Rafikovl (|2010l ) after averaging the Hamiltonian 
of the system (see their Section 3.2), is briefly summarized 
hereinafter since it will help in interpreting the numerical 
results. The authors start their analytical development by 
aligning the constant force along the z-axis and find that 
the system admits 3 integrals of motion 



Ih 



'2n 



(3) 

(4) 



and the semimajor axis, where K(t) = ^(1 - eitf) is the 
dimensionless total angular momentum which is not con- 
served, while Kz is conserved since the force is parallel to 
the z-axis. By solving the Hamilton equations, they find 
that K oscillates between a maximum and minimum value 
given by 

-I 1/2 



i + Ki 



1 + Kl~ Ijj 



Ki 



(5) 



The maximum and minimum values of the angular momen- 
tum translates into minimum and maximum values of ec- 
centricity which depend on the inclination of the ISM fiux 
respect to the initial orbit. The time dependence of K is 
periodic and given by 



Kl + Kit 



+ 



2 2 
The period of the cycle Tstark is 



,(4.i^) (6) 
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(7) 



where a is the orbital semimajor axis, G the gravita- 
tional constant, Mq the star mass and S is the constant 
force applied in the z-axis direction. Tstark is constant in the 
general Stark problem, but it will change with time when we 
deal with dust particles under the action of radiation and 
ISM forces because of the inward drift that reduces the semi- 
major axis a and tends to circularize orbits. As a reference, 
the average value of Tstark for the initial sample of 1 /im 
particles in our model is about 2.3 Myr. This means that 
the dust grains will achieve an eccentricity of about 1 on 
a timescale shorter than Tstark/ 2 ~ l.lMyr. When their 
eccentricity is close to 1 they drift quickly inside. 



Once derived Ti' as a function of t, we can easily com- 
pute e{t), u){t) (from the constant Ih) and i{t) (from Hz)- 
The solution is closed by the equation that gives the node 
longitude as a function of K{t) 



Q.[K) — arctan 



K^ 



■K2){Kl-Ki) 



(8) 



Thanks to the secular integrability of the Stark problem 
we can interpret the evolution shown in Fig[2]and the den- 
sity plots in Fig[l] In the planar case, the maximum value 
of eccentricity cm = 1 and the dust parti cles foll ow the evo- 
lution predicted by the equations of the iBelvaev fc Rafikovl 
(|2010l) model in the {e,uj) plane. However, to properly com- 
pare the analytical solution of the Stark problem with our 
numerical results, we have to take into account the differ- 
ence in the definition of inclination. The rule we adopt in 
this paper is that the inclination is counted from the or- 
bital plane by a counter-clockwise rotation around the x- 
axi s. As a consequence, i = 0° corresponds to i = 90° in 
the lBelvaev fc Rafikovl (|201G| ') model. This introduces a 180° 
difference in the definition of the perihelion longitude be- 
t ween our numerical solutio ns and the analytical formalism 
of lBelv aev fc Rafikovl l|2010l ). In Fig|5]the analytical curves 
for 3 different initial conditions in e and to are shown as 
continuos lines with arrows pointing in the direction of the 
time flow. The particles move from an eccentricity of about 
0.4 and, following the countours of constant Hz and Ih, 
evolve towards e = 1 while the pericenter tends towards the 
asymptotic value of 270° . 

This behaviour explains why the shape of the disk in 
the density plots, for low values of r, appears as a super- 
position of two ellipses. Most particles cluster either from 
below or from above around the fixed peri center v a lue o f 
270° . This align ment was also predic ted by ISchereij 



and observed in iManess et al] (|2009|) . Thanks to the ana- 
lytical solution of IBelvaev fc Rafikovl l|2010i ) it is possible 
to predict the angle between the two elliptical shapes that 
stand out in Fig[T] The analytical curves that enclose, for 
large eccentricities, all the particles, give for e = 1 a value 
of ~ 244.5° and uj ~ 295.5°. Since most of the parti- 
cles cluster close to these two curves, the angular separation 
between the two ellipses of Fig[T] is approximately 51°. A 
shortcut to obtain t he angula r separatio n is t o use directly 
the formula (33) of IBelvaev fc Rafiko^^ (|2010l ') that, in the 
planar case, links eccentricity and perihelion to their con- 
stant initial values through the equation ecosu — A where 
A is the value of ecosu at t=0. By checking Fig|2] we no- 
tice that the curves leading to the outer edges of the peak 
in the (e, to) plane approximately evolve from the following 
points: (0.43,180°) and (0.43,360°). Using these as start- 
ing values to compute A, we get A = ±0.43. As a con- 
sequence, the asymptotic values of u, when e approaches 
1, are ui ~ arcos(0.43) ~ 64.5° + 180° = 244.5° and 
cj - r-cos(-0.43) ~ 295.5° + 180° = 295.5°. The same value 
of separation of ~ 51° is obtained. 

When the eccentricity becomes large, particles begin to 
drift inwards towards the star at increasing speed because of 
the strong PR drag and ISM perturbations on the semimajor 
axis. There is a larger variation in the difference between the 
orbital and ISM velocity over an orbit, and this causes an 
increasingly faster inward drift of grains. This explains why 
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the size of the disk, for small optical depths, is smaller even 
if more asymmetric. 

In our model, when the particles have a semimajor axis 
smaller than 20 AU they are removed from the representa- 
tive sample and a new particle is drawn. This limit is re- 
lated to the large eccentricity of these particles when they 
have drifted inside. With a semimajor axis of 20 AU and 
an eccentricity of 0.9 they are very close to the star. It is 
expected that in this situation the dust particles may ei- 
ther s ublimate or be scattered by inte r nal pl anets. Accord- 
ing to iKobavashi et all (|2009l '): iMukail (liggd ) micron-sized 
dust grains consisting of a silicate core with an icy man- 
tle may undergo active sublimation at temperatures higher 
than 110 K. This temperature is estimated to be present 
on the surface of a 1 /^m grain at about 20 AU from a so- 
lar type star. This scenario, however, does not easily apply 
to dust particles pe rturbed by ISM. They are on highly ec- 
centric orbits while IKobavashi eTal (|2009l ) consider grains 
on circular orbits and perturbed by PR drag only. In our 
model, the drift timescale is very short since, in addition to 
PR drag, a strong contribution to the drag force is given by 
the ISM perturbations and we expect that the particle have 
a rapid infall towards the star when an eccentricity larger 
than 0.9 is achieved. It appears a complex issue to predict 
the amount of sublimation in this context and compare it 
to the drift timescale. For simplicity, we assume that when 
the semimajor axis of a dust particle crosses the lower limit 
of 20 AU a series of mechanisms would make its contribu- 
tion to the outer disk negligible. A test run where this limit 
was lowered to 15 AU did not change the outcome in terms 
of density plots. Scattering by planets appear also a rea- 
sonable outcome when particles move close to the star. The 
presence of a debris disk is a strong indication that planetary 
formation occurred in the system. It is then expected that 
a planetary system orbits the star at some distance from it. 
The potential planets would intercept the dust particles as 
they get close to the star eliminating them from the disk 
population. 

When we consider larger particles, a different balance 
between perturbative changes due to ISM and lifetime is 
reached. In addition, the initial disk is expected to be less 
radially spread since the value of /3 is smaller and the initial 
orbital elements of the grains are closer to those of the par- 
ent bodies. In Fig|3]we show density plots of a dusty disk 
made of 10 /im size particles (/? = 0.03). If the optical depth 
is set to r = 10"'^ the disk does not show any difference from 
the case where ISM is neglected since the ISM perturbations 
do not have enough time to build up a consistent variation 
of eccentricity and perihelion. It appears as a regular ring 
almost overlapping with the parent body ring. The initial 
eccentricities are close to because of the small value of /3 
so that the disk appears circular and the overdense rings ob- 
served in Fig[l] top left plot, are not present. For a smaller 
optical depth, r = 10~^, eccentricities as large as 0.15 are 
pumped up by the ISM perturbations and the disk appears 
slightly eccentric. When finally r = 10~®, the orbits of the 
grains reach a maximum eccentricity of ~ 0.8 and a signifi- 
cant clump develops close to the apocenter of the orbits. In 
this case, the two-ellipse pattern, observed for 1 /xm in Fig[l] 
bottom-right plot, collapses into a single elliptical structure. 
The explanation for this different behaviour resides in the 
lower initial eccentricity with which 10 /^m are produced. 




0.2 0.4 0.6 0.8 1 
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Figure 2. Evolution in the (e, ro) plane of 100 sampled dust 
particles. The orbital elements are computed over an interval of 
time of 5 Myr. The optical depth r is set to 10^® to grant enough 
time to the grains to dynamically evolve. 

This is about 0.03, small compared to the average value of 
about 0.4 for 1 ^m particles. By inspecting Fig|2]we notice 
that the analytical curve of the Stark problem for particles 
passing through low eccentricities gives values of perihelion 
very close to 270° when the eccentricity approaches 1. All 
grains in the 10 /im case evolve from low eccentricity values 
and when their eccentricity is pumped up by the ISM pertur- 
bations, their pericenter is almo st perfectly aligned to 270 . 
Exploiting again formula (33) of iBelvaev fc Rafikovll|2010l ). 
when e = 0.03 , we get an angular separation between the 2 
streams of particles, approaching uj = 270° for large eccen- 
tricities, of only ~ 4°. As a consequence, in the density plot 
the two ellipses merge in a single one. 

The behaviour shown in Fig|3] confirms that the optical 
depth r is the critical parameter in determining the amount 
of asymmetry produced by the ISM flux on a dust disk even 
for larger sizes of the dust particles. When r is large, the 
effects of the ISM flux become negligible and the debris disk 
appears symmetric. 



5 THE INCLINED CASE 

When the ISM fiow is inclined of an angle iisM respect to 
the parent body plane, the analytical model of the Stark 
problem predicts an increase of the grain orbital inclination 
and a reduction of the m aximum value of the eccentricity 
ijBelvaev fc Rafikovllioiol l . The 2-D model, described in the 
previous section, continues to be a good approximation only 
when iisM is lower than 10°. For larger values, the ISM 
perturbations lead to 3-dimensional effects that alters the 
density distribution patterns observed in Fig[l] and, in ad- 
dition, cause significant warping of the disk. However, if the 
optical depth r is of the order of 10"'^, the coUisional dis- 
ruption prevents the building up of significant inclination. 
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Figure 3. 2-D density distribution of 10 /xm size dust particles when a steady state is reached. The top-left plot shows the case with 
optical depth r = 10""^, the top-right plot that with r = 10~*. The last and more perturbed case shown in the bottom plot corresponds 
to an optical depth of r = 10"^. This fi gure has to be compared with Fig[T] 



as it does for the eccentricity, and the warping of the disk is 
negligible. 

This is shown in Fig[3] where the outcome of a simula- 
tion with iisM ~ 60° and r = 10~^ is illustrated as a density 
plot in the (x — z) plane. Only small out-of-plane features 
can be detected in the figure which cannot be interpreted 
neither as warping or clumping due to their limited exten- 
sion along the z-axis (the scale along the z-axis is 1/10 of 
the scale along the x-axis). A large optical depth prevents 
then a disk to develop structures not only in the parent 
body plane but also out-of-plane. This finding is confirmed 



for any value of iisM- In Fig[5]the median inclination of the 
grains in the disk is plotted for different values of iisM with 
r — (top plot). The value of the median of the particle 

inclinations is always smaller than 5° since the Stark cycle is 
quickly interrupted by a collision. This confirms that if the 
disk is optically thick, the ISM does not have time to build 
up observable signatures in the density distribution of de- 
bris disks. In addition, the particles do not have the time to 
drift inside since they are destroyed before PR drag and the 
ISM drag force them to migrate inside. If a gap is present 
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Figure 4. Density distribution in the {x z) plane when r = 

10~^. The quick collisional disruption prevents large out-of-plane 
structures to develop. Notice that the scale along the z-axis is 
1/10 of the scale along the x— axis. 



in the inside regions of the disk this will not be filled up by 
the grain migration. 

Significantly different is the situation when r is as low 
as 10~® (FigO bottom plot): particles in the disk have time 
to move along the Stark cycle. This is not interrupted by 
collisions since the collisional lifetime is longer than Tstark 
and the particles can build up large orbital eccentricities 
and inclinations. In this scenario, the only mechanism able 
to halt the Stark cycling is the fast inward drift due to the 
large eccentricity. The particles migrate inside driven by the 
PR and ISM drag forces ending either onto the star or sub- 
limating or impacting onto a planet. This sink mechanism 
is more effective for smaller values of iisM since the maxi- 
mum eccentricity Emax that can be achieved by the particles 
decreases for increasing iisM, according to the Stark gen- 
eral problem. In Fig[S]we show the evolution of the eccen- 
tricity and inclination for increasing values of iisM- Their 
values are strongly correlated, as predicted by the theory of 
iBelvaev fc Rafikovl (|2010t ) (the theoretical curves are plotted 
as continuous lines in the figures). 

Two interesting features come out from Fig[51 First of 
all the drag forces tend to reduce the eccentricity detaching 
the numerical data from the theoretical curve even if the 
difference is not very marked. In addition, when the max- 
imum eccentricity is lower, for higher iisM, the particles 
move farther within the Stark cycle because their drift time 
is longer. When iisM = 90° the particles complete a full cy- 
cle and some start a new one, even if they do not complete 
more than two cycles before drifting inside 20 AU. 

Additional features related to the dynamics of the dust 
particles under the action of the ISM flow and PR drag can 
be seen in the density plots shown in Fig[7l In the parent 
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Figure 5. Median inclination in a dust disk of 1 l/xm size particles 
for different values of ijSM and r. In the top plot r = 10"'^ while 
in the bottom one it is t = 10~®. 



body plane the density distribution appears very asymmet- 
ric. In the case where iisM = 20° (FigO top plots) it is 
still possible to recognize the two elliptical structures pro- 
duced by the evolution in the (e,a;) plane but their shape is 
blurred due to the inclination distribution. However, when 
iiSM = 60° the double-elliptical structure has fully disap- 
peared and the inclination distribution determines the den- 
sity distribution also in the parent body plane. As it can 
be argued from Fig[6] middle plots, the orbital inclination 
is concentrated around i ~ 50° and, in addition, when the 
inclination is higher, the node longitude is clustered around 
Q. = 7r/2. This leads to an eccentric disk structure that is 
inclined respect to the original orbital plane of the parent 
body of approximately 50°. Finally, the case iisM = 90° 
(Fig[6] bottom plots) is the most symmetric case. The cir- 
cular overdense shape of the disk close to the star is where 
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1. Correlation between the eccentricity and inclination 
(articles in a disk with different values of ijSM- In the 
'SM = 20°, in the middle ijSM = 60° and in the bottom 
0°. The continuous lines are the theoretical predictions 
om the analytical theory of the Stark problem. 
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Figure 8. Eccentricity-inclination distribution for 10 jira size 
dust particles in a disk with iisM = 60°. 



the inclined orbits cross the parent body plane. The dust 
particles have a statistical uniform distribution of the nodes 
for any value of inclination and, as a consequence, a high 
density region is formed in correspondence to the intersec- 
tion between the inclined orbits and the initial parent body 
plane. 

For larger dust grains (10/im) and low optical depth 
(r = 10^®), the dynamical scenario is slightly different. The 
Stark period is of the order of 23 Myr, 10 times longer than 
that of Ifim particles, but the drift rate has not decreased 
by the same amount. As a consequence, the grains evolve 
less into the Stark cycle and migrate inside before reaching 
large values of inclination. This is shown in Fig[8] where the 
eccentricity and inclination of the dust, when the disk is in a 
steady state, are illustrated. The initial eccentricity is close 
to since j3 is small and this leads to a different Stark cycle 
compared to that in Fig[6] middle plots. The particles slowly 
evolve towards larger inclination values but they plunge into 
the inner regions of the disk and are lost before they can de- 
velop significant inclinations. The disk is then expected to be 
less perturbed compared to the case for 1 /im size particles. 
This is confirmed by Fig[9] where the spatial density distri- 
bution is illustrated and can be compared to that shown in 
FigEl (middle plot) for 1 /im size particles. 



6 CONCLUSIONS 

We have shown in this paper, by numerical modeling the 
evolution of debris disks under the effects of solar radiation 
pressure, PR drag and ISM flux, that for typical values of 
optical depth r ~ 10~^ the signatures of the ISM wind on 
grains with radius in the range 1 — 10 fim, just above the cut- 
off size imposed by radiation pressure, are almost negligible. 
The amount of perturbations due to the interaction of dust 
particles with the local ISM neutral atom flow are strongly 
reduced when the grain lifetime is short as in disks with large 
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Figure 7. Density distributions of 1 fim size dust particles in the {x y) and (x z) plane for different values of ijSM- The top 

plots refer to the case iisM = 20°, the middle plots to ijSM = 60° and the lower plots to ijSM = 90°, respectively. 
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Figure 9. Density distributions of 10 /tm size dust particles in the (x y) and {x z) plane for iisM = 60°. Tliese distributions 

can be compared to those for 1 /jm shown in Fig|6] middle plots. 



values of r. In this scenario, the presence of asymmetries 
must be ascribed to different mechanisms like the presence 
of massive bodies within the disk. The neutral ISM fails in 
producing either warping or clumping in such disks. 

Different is the scenario when the optical depth is small. 
We have shown that for r ~ 10"" significant asymmetries 
appear on the density profile of the disk both in the par- 
ent bodies plane and out-of-plane. Interesting is a double- 
elliptical pattern that develops in the parent bodies plane 
when the ISM flow is almost coplanar to it. The disk struc- 
tures are caused by the dynamical evolution of the orbital 
elements of the dust grains. The way in which these orbital 
parameters evolve can be in part predicted on the basis of 
the general Stark model that helps in understanding the 
periodic nature of eccentricity, inclination and the dynam- 
ically related angles perihelion argument an d nodal longi- 
tude jBelvaev fc Rafikovllioiol : iPastoJ I2OIII ). In addition, 
the semimajor axis has a fast inward drift due to the com- 
bination of PR drag and interaction with the ISM. The two 
drag forces, in particular that related to the ISM flow, are 
strong because of the large ecceentricity of the grains. The 
particles quickly migrate towards the star where they can ei- 
ther sublimate or be destroyed by collisions with planets or 
the star itself. This is an additional powerful sink mechanism 
for debris disks with small optical depth. On the contrary, if 
the disk has a large optical depth, coUisional disruption acts 
on a much shorter timescale and it prevents a significant in- 
ward migration. The disk shape in this case would be mostly 
due to the initial value of /3 which sets the initial orbital ele- 
ment distribution of the grains. The debris disk would then 
extend outwards respect to the parent body ring in spite of 
the strong ISM and radiation drag force. This additional dif- 
ferent dynamical behaviour of the dust grains, which again 
depends on the optical depth, sets a nice correlation between 
the value of r and location of the disk respect to the parent 



body ring. Disks with low values of r would expand inward 
respect to the parent body ring down to the sublimation 
region or where planets are orbiting, while disks with large 
values of r would extend mostly outwards. 

Our study suggests that observational data of debris 
disks need care to be interpreted. Potential asymmetries 
identified in the density distribution may be ascribed to 
interactions with the local flux of ISM neutral atoms sur- 
rounding its parent star only if the disk has a low optical 
depth. Otherwise, alternative explanations, like the presence 
of planets, must be investigated. 
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